Bismark软件使用入门

01

甲基化文库建库原理

借助重亚硫酸盐（Bisulfite）转换的方法被认为是金标准，该化学试剂会将未甲基化的C进行转化，胞嘧啶（C）转换为尿嘧啶（U），后续的PCR扩增过程U会被胸腺嘧啶（T）替代，而，甲基化的C则不受影响。在了解如何使用Bismark之前，需要先了解DNA甲基化文库的结构。

通常，甲基化文库构建方法有两种：

• 一种是先片段化，然后末修，加上特殊处理的接头（甲基化的）后，进行CT转化，之后进行扩增，如下图：

该文库为链特异性文库，使用bismark默认参数比对即可。

• 另一种是先CT转换，然后扩增加接头，如下图。

该文库为非链特异性文库，需要使用PBAT文库参数。

通常而言，甲基化建库均为链特异性方式构建：

• 测序read1都是来自于原始的基因组，经历了C->T转化，来自于正链的叫OT，负链叫OB；

• 测序read2都是来自于原始基因组的互补链，经历了G->A转换，称之为CTOT或者CTOB。

02

bismark比对分析

bismark的分析分为三步：

1

基因组索引构建：需要将基因组进行相应的转换，模拟bisulfite处理后基因组，这样才能用于比对。在每个新物种比对之前，都需要进行该处理。bismark内部默认调用bowite软件进行比对。

2

序列比对：输入为下机序列、基因组和比对参数，bismark会输出比对结果和甲基化检测的结果，默认格式为BAM，同时还会有一个统计报告。

3

甲基化位点提取：这一步是可选的，会利用bismark的比对结果来获得甲基化信息。这一步会把甲基化分为不同的类别（CG，CHG和CHH），获得链特异性信息，并且提供过滤参数；也可以对甲基化信息进行调整，或者进行更多的深入分析。

01

基因组索引 

使用bismark_genome_preparation来对基因组进行处理，输入是给定的目录，目录里面是.fa或者.fasta文件。

bismark会生成两个目录，一个是C->T转换的基因组，一个是G->A转换的基因组共两套基因组，之后可以用bowtie2-build来建立索引。

特别注意，索引路径中必需包含一个基因组文件。

运行命令示例如下：

bismark_genome_preparation
    --path_to_bowtie /path/bowtie2/  这里使用bowtie2
    <path_to_genome_folder>   索引路径
     --verbose   输出log信息

索引输出目录结果如下：

Bisulfite_Genome
├── CT_conversion
│   ├── BS_CT.1.bt2
│   ├── BS_CT.2.bt2
│   ├── BS_CT.3.bt2
│   ├── BS_CT.4.bt2
│   ├── BS_CT.rev.1.bt2
│   ├── BS_CT.rev.2.bt2
│   └── genome_mfa.CT_conversion.fa
└── GA_conversion
    ├── BS_GA.1.bt2
    ├── BS_GA.2.bt2
    ├── BS_GA.3.bt2
    ├── BS_GA.4.bt2
    ├── BS_GA.rev.1.bt2
    ├── BS_GA.rev.2.bt2
    └── genome_mfa.GA_conversion.fa

02

比 对 

比对的过程需要提供两个信息

• 输入目录：已经建立好的索引目录，即基因组索引步骤中的索引目录

• 输入文件：测序reads，fastq格式

• 其余参数都是可选，自行调整使用

输出包括：

• 比对的结果，BAM文件

• 报告文件：包含各种统计信息和甲基化统计结果

2.2.1、比对原理

一般情况下，由于建库的时候会产生四种可能性的序列，原始的正链（OT），原始的互补链（OB），原始正链的互补链（CTOT），原始互补链的互补链（CTOB），因此bismark会进行进行四条路线的比对，如下图。

• 左边的C->T，是原始链的转化方式，然后跟两种参考基因组进行比对，得到OT和OB的结果；

• 右边的G->A，是扩增链的转化方式，然后跟两种参考基因组进行比对，得到CTOT和CTOB的结果；

上面示意图提到建库方式是有方向性的，因此bismark只需要进行两次比对即可，而不需要进行4次的比对，如果文库是非链特异性文库，需要加上--non_directional的参数，bismark会进行4次比对，选择最优的结果进行输出。

针对先转化后加接头的文库（PBAT），常见如单细胞甲基化文库，需要添加--pbat参数进行比对。

2.2.2、运行bismark

bismark 
    --bowtie2  使用bowtie2
    -N 0       允许错配数目
    --un       过滤多处匹配reads
    -o         输出目录
    <genome_folder>  索引路径
    -1 read_1.fq.qz  测序read1
    -2 read_2.fq.gz  测序read2

03

输 出 

bam文件：可以利用SeqMonk进行可视化；也可以用于继续的去除dup等

比对输出路径示例如下：

# 样本：A
A_P_R1_bismark_bt2_pe.bam
A_P_R1_bismark_bt2_PE_report.txt
A_P_R2.fq.gz_unmapped_reads_1.fq.gz
A_P_R2.fq.gz_unmapped_reads_2.fq.gz

比对输出路径中包含一个比对后的文件（BAM格式）和一个比对报告文件以及未比对到参考基因组上的reads文件（--un参数）。

若测序数据为双端数据文件中会有PE或pe标签，若为单端数据则会有SE或se字符加以区别。

针对链特异性文库，若进行单端比对read1使用bismark单端参数即可，read2需要添加--pbat参数；而非链特异性文库--pbat参数的使用则相反。

通常，若某物种比对率较低，可先进行双端比对同时输出未必对的reads，随后使用未比对的read1和read2分别进行单端比对并合并输出的BAM文件，作为下游的输入。

比对报告文件示例如下：

Bismark report for: A_P_R1.fq.gz and A_P_R2.fq.gz
Bismark was run with Bowtie 2 against the bisulfite genome of /Index/ with the specified options: -q --score-min L,0,-0.2 --reorder --ignore-quals --no-mixed --no-discordant --dovetail --maxins 1000
Option '--directional' specified (default mode): alignments to complementary strands (CTOT, CTOB) were ignored (i.e. not performed)

Final Alignment report
======================
Sequence pairs analysed in total:  43930319
Number of paired-end alignments with a unique best hit: 26314484
Mapping efficiency:  60.0%
Sequence pairs with no alignments under any condition:  16147263
Sequence pairs did not map uniquely:  1468572
Sequence pairs which were discarded because genomic sequence could not be extracted:    2

Number of sequence pairs with unique best (first) alignment came from the bowtie output:
CT/GA/CT:   13199639   ((converted) top strand)
GA/CT/CT:   0  (complementary to (converted) top strand)
GA/CT/GA:   0  (complementary to (converted) bottom strand)
CT/GA/GA:   13114843  ((converted) bottom strand)

Number of alignments to (merely theoretical) complementary strands being rejected in total:     0

Final Cytosine Methylation Report
=================================
Total number of C's analysed: 1191141867

Total methylated C's in CpG context: 33527315
Total methylated C's in CHG context: 17696056
Total methylated C's in CHH context: 68376839
Total methylated C's in Unknown context: 97299

Total unmethylated C's in CpG context: 20813883
Total unmethylated C's in CHG context: 239796745
Total unmethylated C's in CHH context: 810931029
Total unmethylated C's in Unknown context:  301701

C methylated in CpG context: 61.7%
C methylated in CHG context: 6.9%
C methylated in CHH context: 7.8%
C methylated in unknown context (CN or CHN): 24.4%

Bismark completed in 0d 7h 30m 40s

报告展示了唯一比对上的read对数目等信息（Final Alignment report部分）以及有reads支持发生甲基化的C的数目统计（Final Cytosine Methylation Report部分），注意这里并非真正意义上的甲基化C，仍需要下游的分析。

04

去除duplicated序列

由于PCR的扩增，可能会存在比对到基因组相同位置的序列，重复测序的序列需要被去除。如果是RRBS文库，建议不要进行这一步骤。

对于单选测序序列，去dup的时候，会根据染色体、起始坐标和链的方向来确定是否是duplication；

对于双端测序序列，Read1的染色体、起始坐标，Read2的起始坐标会用于去除duplication。需要注意的是，去dup的时候需要注意read1和read2是相邻的，而不是按照坐标排序的。

去除dup的参数示例如下：

deduplicate_bismark 
    -p A_P_R1_bismark_bt2_pe.bam   bismark比对输出的bam文件 
    --bam                          bam格式输出
    --samtools_path <samtools_dir> SAMtools软件路径
    --output_dir <out_dir>         输出路径

输出路径会生成一个去除dup后的bam文件和一个去除dup的统计文件。

03

提取甲基化位点

使用bismark_methylation_extractor来提取每个C位点的甲基化的情况，结果文件可以导入到SeqMonk中。

可以使用--bedGraph，来产生bedgraph和coverage的文件；--count可以产生count文件，输出每个CpG位点的深度文件。

bismark_methylation_extractor
    --pair-end         指定双端数据
    --comprehensive    输出CHG CHH CpG的甲基化信息
    --no-overlap       去除reads重叠区域的bias
    --bedGraph         输出bedGraph文件
    --counts           每个C上甲基化reads和非甲基化reads的数目
    --report           一个甲基化summay
    --cytosine_report  输出全基因组所有CpG
    --genome_folder    <path_to_reference_genome>
    input.bam          输入文件
    -o >output_dir>    输出路径

bismark对不同状态的胞嘧啶进行了约定并记录到了BAM文件和中间文件中：

X 代表CHG中甲基化的C
x  代笔CHG中非甲基化的C
H 代表CHH中甲基化的C
h  代表CHH中非甲基化的C
Z  代表CpG中甲基化的C
z  代表CpG中非甲基化的C
U 代表其他情况的甲基化C(CN或者CHN)
u  代表其他情况的非甲基化C (CN或者CHN)


bismark_methylation_extractor的一个非常重要的输出文件为全基因组胞嘧啶报告文件（*CX_report.txt），该文件是一个非常核心的记录胞嘧啶深度信息的文件，可进行甲基化水平的计算、区域甲基化信号文件的转化、差异甲基化水平的计算等。

文件示例如下：

chr1    14716   +       1       6       CG      CGT
chr1    14717   -       10      1       CG      CGA
chr1    14741   +       2       1       CG      CGG
chr1    14742   -       0       1       CG      CGC

各列信息解释如下：

第1列 : 染色体编号
第2列 : 染色体起始位置
第3列 : 正负链信息
第4列 : 甲基化碱基数目
第5列 : 非甲基化碱基数目
第6列 : CG/CHG/CHH
第7列 : C背景序列

基于该文件可进行单个胞嘧啶位点的甲基化水平（C/C+T）计算，比如第一行：甲基化是水平为1/(1+6)；针对指定区域的甲基化水平计算则取区域内的胞嘧啶的总C/(总C+总T)即可。

04

bismark的优缺点

4.1、优点

（1）比对和mC Calling比较准确，受众广泛，使用率很高

（2）软件维护性好，更新快，bug修复及时

（3）软件成体系，集成了比对、去重、mC Calling分析

4.2、缺点

（1）比对和mC Calling速度慢，真的很慢，可以通过拆分输入数据为小份作为输出和增大计算资源解决。

（2）bismak相匹配的直接进行下游甲基化水平、图形可视化等分析的套件较少，需要自行计算或进行格式转化作为第三方软件的输入。

（3）说明文档不够系统和详细，在”外行“看来理解从建库到mC Calling有一定的挑战。

bismark的结果只是研究甲基化信息的第一步也是关键的一步。未来，期望bismark在速度上，和在对甲基化水平层面分析的支持上多一下改进。

05

参考资料

1.https://github.com/FelixKrueger/Bismark

2.https://github.com/FelixKrueger/Bismark/tree/master/Docs

3.https://www.epigenesys.eu/images/stories/protocols/pdf/20120720103700_p57.pdf

作者：马可菠萝 童蒙

编辑：angelica

往期精彩回顾